//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "InterfaceDiffusion.h"

registerMooseObject("snakeApp", InterfaceDiffusion);

template <>
InputParameters
validParams<InterfaceDiffusion>()
{
  InputParameters params = validParams<InterfaceKernel>();
  params.addRequiredParam<Real>("sigma", "sigma");
  params.addRequiredParam<Real>("epsilon", "epsilon");
  params.addRequiredParam<Real>("masterkappa", "The coefficient of diffusion");
  params.addRequiredParam<Real>("pairedkappa", "The coefficient of diffusion");
  return params;
}

InterfaceDiffusion::InterfaceDiffusion(const InputParameters & parameters)
  : InterfaceKernel(parameters),
     _epsilon(getParam<Real>("epsilon")),
    _sigma(getParam<Real>("sigma")),
    _ownkappa(getParam<Real>("masterkappa")),
    _neikappa(getParam<Real>("pairedkappa"))
{
}

Real
InterfaceDiffusion::computeQpResidual(Moose::DGResidualType type)
{
  Real r = 0;
  
  const unsigned int elem_b_order = _var.order();
  double h_elem =
      _current_elem->volume() / _current_side_elem->volume() * 1. / Utility::pow<2>(elem_b_order);

  h_elem /= _ownkappa;
  switch (type)
  {
    case Moose::Element:
      r -= 0.5 * (_ownkappa * _grad_u[_qp] * _normals[_qp] +
                  _neikappa * _grad_neighbor_value[_qp] * _normals[_qp]) *
           _test[_i][_qp];
      r += _epsilon * 0.5 * (_u[_qp] - _neighbor_value[_qp]) * _ownkappa * _grad_test[_i][_qp] *
           _normals[_qp];
      r += _sigma / h_elem * (_u[_qp] - _neighbor_value[_qp]) * _test[_i][_qp];
      break;

    case Moose::Neighbor:
      r += 0.5 * (_ownkappa * _grad_u[_qp] * _normals[_qp] +
                  _neikappa * _grad_neighbor_value[_qp] * _normals[_qp]) *
           _test_neighbor[_i][_qp];
      r += _epsilon * 0.5 * (_u[_qp] - _neighbor_value[_qp]) * _ownkappa *
           _grad_test_neighbor[_i][_qp] * _normals[_qp];
      r -= _sigma / h_elem * (_u[_qp] - _neighbor_value[_qp]) * _test_neighbor[_i][_qp];
      break;
  }

  return r;
}

Real
InterfaceDiffusion::computeQpJacobian(Moose::DGJacobianType type)
{
  Real r = 0;
  
  const unsigned int elem_b_order = _var.order();
  double h_elem =
      _current_elem->volume() / _current_side_elem->volume() * 1. / Utility::pow<2>(elem_b_order);

  switch (type)
  {
    case Moose::ElementElement:
      r -= 0.5 * _ownkappa * _grad_phi[_j][_qp] * _normals[_qp] * _test[_i][_qp];
      r += _epsilon * 0.5 * _phi[_j][_qp] * _ownkappa * _grad_test[_i][_qp] * _normals[_qp];
      r += _ownkappa*_sigma / h_elem * _phi[_j][_qp] * _test[_i][_qp];
      break;

    case Moose::ElementNeighbor:
      r -= 0.5 * _neikappa * _grad_phi_neighbor[_j][_qp] * _normals[_qp] * _test[_i][_qp];
      r += _epsilon * 0.5 * -_phi_neighbor[_j][_qp] * _ownkappa * _grad_test[_i][_qp] *
           _normals[_qp];
      r += _ownkappa*_sigma / h_elem * -_phi_neighbor[_j][_qp] * _test[_i][_qp];
      break;

    case Moose::NeighborElement:
      r += 0.5 * _ownkappa * _grad_phi[_j][_qp] * _normals[_qp] * _test_neighbor[_i][_qp];
      r += _epsilon * 0.5 * _phi[_j][_qp] * _neikappa * _grad_test_neighbor[_i][_qp] *
           _normals[_qp];
      r -= _neikappa*_sigma / h_elem * _phi[_j][_qp] * _test_neighbor[_i][_qp];
      break;

    case Moose::NeighborNeighbor:
      r += 0.5 * _neikappa * _grad_phi_neighbor[_j][_qp] * _normals[_qp] *
           _test_neighbor[_i][_qp];
      r += _epsilon * 0.5 * -_phi_neighbor[_j][_qp] * _neikappa *
           _grad_test_neighbor[_i][_qp] * _normals[_qp];
      r -= _neikappa*_sigma / h_elem * -_phi_neighbor[_j][_qp] * _test_neighbor[_i][_qp];
      break;
  }

  return r;
}
